Extracción de reflectancia

La reflectancia de superficie (R_rs) se extrae usando las coordenadas geográficas de los sitios de muestreos y el recorte del producto satelital, luego de ser procesado por ACOLITE.

Recorte del área de estudio

ACOLITE genera un archivo ráster por cada banda espectral de S2-MSI, ya recortados al área de interés del Proyecto, mostrada a continuación.

In [1]:

v <- terra::vect("../vector/roi_acolite.geojson") |>
  terra::project("EPSG:4326")

ce <- terra::centroids(v) |> 
  terra::geom() |> 
  tibble::as_tibble()

leaflet::leaflet() |>
  leaflet::addProviderTiles(
    provider = "Esri.WorldImagery"
  ) |> 
  leaflet::addPolygons(
    data = v,
    color = "#B86092",
    stroke = TRUE,
    fillColor = "transparent",
    opacity = 1,
    weight = 5
  ) |>
  leaflet.extras::addFullscreenControl(
    position = "bottomright"
  ) |>
  leaflet.extras::addResetMapButton() |> 
  leaflet::setView(zoom = 10, lat = ce$y, lng = ce$x)

Vector utilizado para el recorte del producto.

Sitios de muestreo

Las coordenadas de los sitios de muestro se extraen de los datos registrados durante las campañas de campo.

In [2]:

coord_sitios <- readxl::read_xlsx(
    path = "datos/datos_gistaq.xlsx",
    sheet = 1,
    .name_repair = "unique_quiet"
  ) |>
  select(fechas = 1, longitud = 4, latitud = 5) |>
  fill(fechas) |>
  mutate(fechas = ymd(fechas)) |>
  dplyr::filter(fechas == ymd(fecha_i)) |>
  drop_na() |>
  select(-fechas) |>
  mutate(punto = row_number(), .before = 1)

coord_sitios_v <- vect(
  coord_sitios,
  geom = c("longitud", "latitud"),
  crs = "EPSG:4326",
  keepgeom = FALSE
) |>
  terra::project("EPSG:32721")

1

Indico la ubicación del archivo Excel.

2

Selecciono la fecha de interés.

3

Cambio las coordenadas para ser coincidente con el producto satelital.

Extracción de los valores de píxel

Leo los rásters almacenados, únicamente de la fecha de interés. Luego, modifico los valores de píxel a promedios de una ventana de 3x3. Con las coordenadas de muestreo, se extraen las reflectancias corregidas por ACOLITE.

In [3]:

r <- list.files(
    path = "recorte_acolite/",
    pattern = as.character(fecha_i),
    full.names = TRUE
  ) |>
    rast()

r_3x3 <- terra::focal(
  r,
  w = 3,
  fun = mean,
  na.rm = TRUE
)

reflect <- terra::extract(r_3x3, coord_sitios_v) |>
  as_tibble() |>
  rename(punto = ID) |>
  pivot_longer(
    cols = -punto,
    names_to = "banda",
    values_to = "reflect"
  ) |>
  mutate(fecha = ymd(fecha_i), .before = 1) |>
  inner_join(coord_sitios, by = join_by(punto))

1

Carpeta donde se almacenan los rásters generados por ACOLITE.

2

Filtro los ráster para procesar solamente los de la fecha de interés.

3

Genero una ventana móvil de 3x3, calculando el promedio de los píxeles.

4

Realizo la extracción de los valores.

5

Combino las reflectancias con las coordenadas de los sitios de muestreo.

Almaceno las reflectancias de superficie.

In [4]:

base_de_datos <- read_csv(datos, show_col_types = FALSE)

bind_rows(base_de_datos, reflect) |>
  arrange(fecha, punto) |>
  write_csv(datos)

1

Leo la base de datos existente.

2

Combino los nuevos datos extraídos con los anteriores.